#**********************************************************************************
# ccamTables.r
# This file contains the code for writing main body tables to disk.
# This file assumes that an object called 'opList' exists and is a valid opList as
#  described in loadScenarios.r
#
# Author            : Chris Grandin
# Development Date  : January 2012 - September 2012
#
#**********************************************************************************

mb.table.pop.est <- function(weightFactor=1,ci=c(0.05,0.95),roundDec=2,writeCSV=T,writeTEX=F,silent=F){
  # makes two tables, mb.table.pop.est and mb.table.pop.est.ci, the
  # summary tables for the population estimates from the MCMC run and the
  # credible interval table for some of those estimates.
  
  years <- A$yrs
  recruitYears <- ryr  # Age-1 recruits different years than the rest

  biomass <- A$mc.sbt
  age3PlusBiomass <- A$mc.bt3

  depletion <- A$mc.sbdepletion
  age1Recruits <- A$mc.rt
  f40spr <- A$mc.sprstatus_f40

	catch <- A$ct/weightFactor
  expFrac <- catch/age3PlusBiomass
  
  pop.est.table <- as.data.frame(matrix(nrow=length(years),ncol=6))
  pop.est.ci.table <- as.data.frame(matrix(nrow=length(years),ncol=6))
  colnames(pop.est.table) <- c("Year","SpawningBiomass","Depletion","Age1Recruits","F40SPR","ExploitationFraction")
  colnames(pop.est.ci.table) <- c("Year","SpawningBiomass","Depletion","Age1Recruits","F40SPR","ExploitationFraction")

  numNAYearsAge1Recruits <- 0
  for(year in 1:length(years)){
    biomassM <- median(biomass[,year])
    age3PlusBiomassM <- median(age3PlusBiomass[,year])

    biomassCI <- quantile(biomass[,year],ci)
    age3PlusBiomassCI <- quantile(age3PlusBiomass[,year],ci)

    depletionM <- median(depletion[,year])
    depletionCI <- quantile(depletion[,year],ci)

    expFracM <- median(expFrac[,year])
    expFracCI <- quantile(expFrac[,year],ci)
    
    if(years[year] %in% ryr){
      age1RecruitsM <- median(age1Recruits[,year-numNAYearsAge1Recruits])
      age1RecruitsCI <- quantile(age1Recruits[,year-numNAYearsAge1Recruits],ci)
    }else{
      numNAYearsAge1Recruits <- numNAYearsAge1Recruits + 1      
      age1RecruitsM <- NA
      age1RecruitsCI <- NA
    }
    if(ncol(f40spr)<length(years) & year==length(years)){
      f40sprM <- NA
      f40sprCI <- NA
    }else{
      f40sprM <- median(f40spr[,year])
      f40sprCI <- quantile(f40spr[,year],ci)
    }

    printMask <- paste("%02.",roundDec,"f",sep="")
    perPrintMask <- "%3.0f"
    pop.est.table[year,] <- list(years[year],
                                 sprintf(printMask,round(biomassM,roundDec)),
                                 paste(sprintf(perPrintMask,round(depletionM*100)),"%",sep=""),
                                 sprintf(printMask,round(age1RecruitsM,roundDec)),
                                 sprintf(printMask,round(f40sprM,roundDec)),
                                 sprintf(printMask,round(expFracM,roundDec)))

    pop.est.ci.table[year,] <- list(years[year],
                                    paste(sprintf(printMask,round(biomassCI[1],2)),"-",sprintf(printMask,round(biomassCI[2],2)),sep=""),
                                    paste(sprintf(perPrintMask,round(depletionCI[1]*100)),"%-",sprintf(perPrintMask,round(depletionCI[2]*100)),"%",sep=""),
                                    paste(sprintf(printMask,round(age1RecruitsCI[1],2)),"-",sprintf(printMask,round(age1RecruitsCI[2],2)),sep=""),
                                    paste(sprintf(printMask,round(f40sprCI[1],2)),"-",sprintf(printMask,round(f40sprCI[2],2)),sep=""),
                                    paste(sprintf(printMask,round(expFracCI[1],2)),"-",sprintf(printMask,round(expFracCI[2],2)),sep=""))
  }
  if(writeCSV){
    fn <- paste(tabDir,"mb.table.pop.est.csv",sep="")
    write.csv(pop.est.table,file=fn,row.names=F)
    fn <- paste(tabDir,"mb.table.pop.ci.est.csv",sep="")
    write.csv(pop.est.ci.table,file=fn,row.names=F)
  }
  if(writeTEX){
  }
  
}

mb.table.age.est <- function(sAge=1,nAge=15,roundDec=2,formatOut="%1.2f",writeCSV=T,writeTEX=F){
  nt <- A$nt
  colnames(nt) <- sAge:nAge
  roundN <- round(nt,roundDec)
  numAtAge <- as.data.frame(matrix(nrow=nrow(roundN),ncol=ncol(roundN)))
  
  for(row in 1:nrow(roundN)){
    numAtAge[row,] <- sprintf(formatOut,roundN[row,])
  }
  numAtAge <- cbind(A$yr,numAtAge)
  colnames(numAtAge)=c("Year",sAge:nAge)
  
  if(writeCSV){
    fn <- paste(tabDir,"mb.table.age.est.csv",sep="")
    write.csv(numAtAge,file=fn,row.names=F)    
  }
  if(writeTEX){
  }
}

mb.table.mle.vs.post <- function(roundDec=3,formatOut="%1.2f",writeCSV=T,writeTEX=F){
  mle <- A
  mcmc <- A$mc
  
  # Negative log likelihood chunk - 4 items
  numRows <- 5
  # Parameters chunk - 5 items
  numRows <- numRows + 6
  # Reference points chunk 4 items
  numRows <- numRows + 4

  d <- as.data.frame(matrix(nrow=numRows,ncol=3))
  colnames(d) <- c("X","MLE","Posterior median")
  d[,1] <- c("Survey index",
             "Commercial ages",
             "Survey ages",
             "Parameter priors",
             "",
             "R0 (billions)",
             "Steepness (h)",
             "Natural mortality (M)",
             "Acoustic catchability (Q)",
             "Additional acoustic survey SD",
             "",
             "2008 Recruitment deviation",
             "SB0 (million mt)",
             "2011 Depletion",
             "2010 SPR Ratio")

  # Add MLE values
  d[1,2] <- sprintf(formatOut,round(mle$nloglike[1],roundDec))
  d[2,2] <- sprintf(formatOut,round(mle$nloglike[2],roundDec))
  d[3,2] <- sprintf(formatOut,round(mle$nloglike[3],roundDec))
#  d[4,2] <- sprintf(formatOut,round(mle$nloglike_p,roundDec))

  d[5,2] <- ""
  
  d[6,2] <- sprintf(formatOut,round(mle$Ro,roundDec))
  d[7,2] <- sprintf(formatOut,round(mle$steepness,roundDec))
  d[8,2] <- sprintf(formatOut,round(mle$M,roundDec))
  d[9,2] <- sprintf(formatOut,round(mle$q,roundDec))
  d[10,2] <- sprintf(formatOut,round(mle$rho/mle$varphi,roundDec))

  d[11,2] <- ""

  d[12,2] <- sprintf(formatOut,round(mle$wt[length(mle$wt)-1],roundDec))  # 2009 recruitment anomaly in 2011 assessment year
  d[13,2] <- sprintf(formatOut,round(mle$sbo,roundDec))
  d[14,2] <- sprintf(formatOut,round(mle$sbt[length(mle$sbt)]/mle$sbo,roundDec))
  d[15,2] <- sprintf(formatOut,round(mle$spr[length(mle$spr)],roundDec))

  # Add MCMC median posterior values - no negative log likelihoods.
  d[1,3] <- sprintf(formatOut,round(median(mcmc$Ro),roundDec))
  d[2,3] <- sprintf(formatOut,round(median(mcmc$h),roundDec))
  d[3,3] <- sprintf(formatOut,round(median(mcmc$m),roundDec))
  d[4,3] <- sprintf(formatOut,round(median(mcmc$q),roundDec))

  d[5,3] <- ""
  
  d[6,3] <- sprintf(formatOut,round(median(mcmc$Ro),roundDec))
  d[7,3] <- sprintf(formatOut,round(median(mcmc$h),roundDec))
  d[8,3] <- sprintf(formatOut,round(median(mcmc$m),roundDec))
  d[9,3] <- sprintf(formatOut,round(median(mcmc$q),roundDec))
  #d[10,3] <- sprintf(formatOut,round(median(mcmc$),roundDec))
  
  d[11,3] <- ""

  #d[12,3] <- sprintf(formatOut,round(median(),roundDec))
  d[13,3] <- sprintf(formatOut,round(median(mcmc$SBo),roundDec))
  d[14,3] <- sprintf(formatOut,round(median(mcmc$depletion),roundDec))
  #d[15,3] <- sprintf(formatOut,round(median(mcmc$),roundDec))

  if(writeCSV){    
    fn <- paste(tabDir,"mb.table.mle.vs.post.csv",sep="")
    write.csv(d,file=fn,row.names=F)    
  }
  if(writeTEX){
  }
}

mb.table.all.param.est <- function(sAge=2,nAge=15,roundDec=2,formatOut="%1.2f",writeCSV=T,writeTEX=F){
  # lastYear is the year corresponding to the last value in the A$wt (recruitment resids)
  mle <- A
  mcmc <- A$mc

  numParams <- 2
  d <- as.data.frame(matrix(nrow=numParams,ncol=3))
  colnames(d) <- c("Parameter","MLE","MCMC median")
  row <- 1

  d[row,1] <- "MSY"
  d[row,2] <- sprintf(formatOut,round(mle$MSY,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$msy),roundDec))
  row <- row + 1
   
  d[row,1] <- "FMSY"
  d[row,2] <- sprintf(formatOut,round(mle$Fmsy,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$fmsy),roundDec))
  row <- row + 1
 
  d[row,1] <- "M"
  d[row,2] <- sprintf(formatOut,round(mle$M,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$m),roundDec))
  row <- row + 1

  d[row,1] <- "Age at 50% 1st harvest"
  d[row,2] <- sprintf(formatOut,round(mle$ah,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$ahat),roundDec))
  row <- row + 1

  d[row,1] <- "SD fishery selectivity"
  d[row,2] <- sprintf(formatOut,round(mle$ghat,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$ghat),roundDec))
  row <- row + 1

  d[row,1] <- "Observation error"
  d[row,2] <- sprintf(formatOut,round(mle$rho,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$rho),roundDec))
  row <- row + 1

  d[row,1] <- "Inverse total variance"
  d[row,2] <- sprintf(formatOut,round(mle$varphi,roundDec))
  d[row,3] <- sprintf(formatOut,round(median(mcmc$varphi),roundDec))
  row <- row + 1

  years <- A$ryrs

  row <- 8
  tmpMeds <- apply(A$mc.recDevPosts,2,median)
  for(year in 1:length(years)){
    d[row,1] <- paste("Mean recruitment #",year," - ",years[year],sep="")
    d[row,2] <- sprintf(formatOut,round(mle$wt[year],roundDec))
    d[row,3] <- sprintf(formatOut,round(tmpMeds[year],roundDec))
    row <- row + 1
  }
   
  if(writeCSV){    
    fn <- paste(tabDir,"mb.table.est.parameters.csv",sep="")
    write.csv(d,file=fn,row.names=F)    
  }
  if(writeTEX){
  }
}

table.29 <- function(){
	#This is the final summary table for the main body
	#Caption in 2011 doc was: Select parameters, derived quantities and reference point estimates
  #for TINSS sensitivity analyses to the prior for M. Note that recruits
  # are age 1 and not directly comparable with SS.
	lbl <- c("Unfished age-1 recruits (billions)",
           "Steepness (h)",
           "Natural mortality (M)",
           "Acoustic catchability (q)",
           "2008 Recruitment",
           "SB0 (million mt)",
           "2012 Depletion",
           "2011 Fishing intensity (1-SPR/1-SPR40%)",
           #
           "Female spawning biomass (SB40% million mt)",
           "SPRSB40%",
           "Exploitation fraction (ct/Bt3) resulting in SB40%", 
           "Yield at SB40% (million mt)",
           "Female spawning biomass (SBF40% million mt)",
           "SPRMSY-proxy",
           "Exploitation fraction corresponding to SPR ",
           "Yield at SBF40% (million mt)",
           "Female spawning biomass (SBMSY million mt)",
           "SPRMSY",
           "Exploitation fraction corresponding to SPRMSY",
           "MSY (million mt)")

  mcs <- A$mcRefPoints
  post.mcs <- as.data.frame(window(mcmc(mcs),start=Burn,thin=Thin))
  mcsci <- t(apply(post.mcs,2,quantile,probs=c(0.025,0.5,0.975)))
  meds <- mcsci[,2]

  ti <- as.numeric(meds[2])  # R0
  ti <- c(ti,median(A$mc$h[(Burn+1):length(A$mc$h)])) # steepness
  
  ti <- c(ti,median(exp(A$mc$log.m[(Burn+1):length(A$mc$log.m)]))) # natural mortality

  ti <- c(ti,median(A$mc$q[(Burn+1):length(A$mc$h)])) # catchability
  
  rt <- A$mc.rt[,43]
  rt <- rt[(Burn+1):length(rt)]
  ti <- c(ti,median(rt)) # 2009 Age-1 recruits

  ti <- c(ti,median(A$mc$bo[(Burn+1):length(A$mc$bo)]))  # unfished biomass

  depletion <- A$mc.sbdepletion[,ncol(A$mc.sbdepletion)]
  depletion <- depletion[(Burn+1):length(depletion)]
  ti <- c(ti,median(depletion)) # last year depletion
  sprstatus <- A$mc.sprstatus_f40[,ncol(A$mc.sprstatus_f40)]
  sprstatus <- sprstatus[(Burn+1):length(sprstatus)]
  ti <- c(ti,median(sprstatus))  # last year spr40 status

  for(i in c(3:6,8:15)){
    ti <- c(ti,mcsci[i,2])
  }

  ti <- cbind(lbl,ti)
  colnames(ti) <- c("Quantity", "Median")
  fn <- paste(tabDir,"table.13.csv",sep="")  # table 13 in 2012 assessment
  write.csv(ti,file=fn, row.names=F)
}

table.1	<- function(){
  fthin	<-	function(samples){
    window(mcmc(samples),thin=Thin)
  }

	sb <- fthin(A$mc.sbt)
	ft <- fthin(A$mc.ft)
	sb.mci <- t(apply(sb,2,quantile,probs=c(0.5,0.05,0.95))) # spawning stock biomass
	dt.mci <- t(apply(sb/sb[,1],2,quantile,probs=c(0.5,0.05,0.95))) #depletion
	ft.mci <- t(apply(ft,2,quantile,probs=c(0.5,0.05,0.95))) #fishing mortality (not shown in assessment Table a)
	ft.mci <- rbind(ft.mci,rep(NA,3))
	
	t.1	<- cbind(yrs,sb.mci[1:nyrs,],dt.mci[1:nyrs,],ft.mci)
	t.1	<- apply(t.1,2,round,3)
	filename <- paste(tabDir,"table1.tex",sep="")
	filenamecsv <- paste(tabDir,"table1.csv",sep="")
  
	cap <- "Median estimate and 5\\% and 95\\% credible intervals for the female spawning
				stock biomass (million mt), spawning stock depletion, and fishing mortality
				rates in 1966 and recent years.  These estimates are based on sampling the
				joint posterior distribution using MCMC, chain length 2,000,000 with 
				systematic samples drawn every 200 iterations."
	cgrp <- c(" ","Spawning stock biomass","Depletion","Fishing Mortality")
	ncgrp	<- c(1,3,3,3)
	colnames(t.1)	<-c("Year",rep(c("median","5\\%","95\\%"),3))
	latex(rbind(t.1[1,],tail(t.1)),file=filename,rowname=NULL,caption=cap,cgroup=cgrp,n.cgroup=ncgrp,label="TinSS.T1")
  write.csv(t.1,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.2	<- function(){   
	P <- A$proj #na.exclude(read.table("tinss.proj",header=T,sep="\t",fill=T))
  # PROBLEM - no .proj file for ccam...
	x <- 0:24
	y <- seq(1,dim(P)[1],by=25*thin)
	ix <- as.vector(sapply(y,"+",x))
	print(dim(P))
	P <- P[ix,]
	print(dim(P))
	with(P,{
    fn.glm <- function(xx,yy){	#xx is a ratio 0 = good, 1= bad
      yy[yy<=1] <- 0
      yy[yy!=0] <- 1	
      fit <- glm(yy~xx,family=binomial(logit))
      return(fit)
    }
  })
	risk <- seq(0.05,0.95,by=0.05)
	log.risk <- log(risk/(1-risk))
	
	fn.tac <- function(fglm){
    x<-round((log.risk-coef(fglm)[1])/coef(fglm)[2],3); x[x<0]=0; return(x)
  }
  
  tac1 <- fn.tac(dt1)
  tac2 <- fn.tac(dt2)
  tac3 <- fn.tac(dt3)
  tac4 <- fn.tac(dt4)
  tac5 <- fn.tac(dt5)
  t.2	<- cbind(risk,tac1,tac2,tac3,tac4,tac5)
  colnames(t.2) <- c("Risk level","$F_{2010}\\geq $\\fmsy","$SB_{2011}\\leq SB_{2010}$","$SB_{2011} \\leq SB_{MSY}$","$SB_{2011} \\leq SB_{40}$","$SB_{2011} \\leq SB_{25}$")
  filename <- paste(tabDir,"T2.tex",sep="")
  filenameb <- paste(tabDir,"T2b.tex",sep="")
  filenamec <- paste(tabDir,"T2c.tex",sep="")
  cap	<-"Decision table for catch advice. The risk level represents the probability
			   of exceeding a specified management target for a given ABC option. The interpretation of this table is as follows;
			   if the management goal is not to exceed the target fishing mortality rate of \\fmsy in 2009 with a 0.25 probability, then the ABC option
			   should be set at 0.067 million mt or less.  If the management target is prevent further decline in spawning stock biomass  with a 0.5 probability
			   then the ABC should be set at 0.111 million mt or less."
  latex(t.2,file=filename,rowname=NULL,caption=cap,label="TinSS.T2")
  latex(t.2,file=filenameb,rowname=NULL,caption=cap,label="TinSS.T2b") 
  latex(t.2[c(5, 10, 15), ],file=filenamec,rowname=NULL,caption=
        "Decision table for catch advice. The risk level represents the probability of exceeding a specified management target.",label="TinSS.T2c")
  print(head(t.2[, 1:3], 10))
 
}

table.3	<-	function(){
	t.3 <- cbind(yrs,A$bt,A$sbt,A$sbt/A$sbo,c(A$ct,NA)/1e6,c(A$ft,NA),A$bt2,A$bt3,c(A$ct/1e6,NA)/A$bt2,c(A$ct/1e6,NA)/A$bt3)
	colnames(t.3) <- c("Year","$B_t$","$SB_t$","$SB_t/SB_{0}$","$C_t$","$F_t$","$B_{t,2+}$","$B_{t,3+}$","$C_t/B_{t,2+}$","$C_t/B_{t,3+}$")
	filename <- paste(tabDir,"Table3.tex",sep="")
	filenamecsv <- paste(tabDir,"Table3.csv",sep="")
	cap <- "Maximum likelihood estimates of vulnerable biomass ($B_t$), spawning biomass ($SB_t$) and depletion,
			 landings ($C_t$ millions mt), instantaneous fishing mortality rates ($F_t$), 2+ and 3+ biomass
			($B_{t,2+}, B_{t,3+}$), and total catch over 2+ and 3+ biomass ($C_t/B_{t,2+}$, $C_t/B_{t,3+}$), from 1966 
			to the begining of \\nnyr."
	latex(round(t.3,2),file=filename,rowname=NULL,caption=cap,label="TinsSS.T3",longtable=T,lines.page=100)
  write.csv(t.3,file=filenamecsv, row.names=F)
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.4 <- function(){	#Table with catch, mean age and survey abundance index
	#may want to add mean age of the survey comps.
	x <- data.frame(Year=A$yr,Ct=A$ct)
	ii <- A$yr%in%A$ayr
	mu.age <- rep(NA,length(A$yr))
	mu.age[ii] <- 2:15%*%t(A$Adata[,1:14])
	x <- data.frame(x,muage=round(mu.age,2))
	
	ii <- A$yr%in%A$iyr
	yt <- rep(NA,length(A$yr))
	yt[ii] <- A$yt
	x <- data.frame(x,It=yt)
	x <- rbind(x,rep(NA,length=4))
	colnames(x) <- c("Year", "$C_t$", "$\\bar{a}$", "$I_t$")
	Table3 <- cbind(x[1:22,],x[23:44,])
	#Table3=x
	cap <- "Combined historical landings (mt) for the U.S. and Can. fisheries, 
	        mean age of the catch, and survey abundance indices (millions mt) from 
	        the acoustic-trawl survey."
	filename <- paste(tabDir,"table4.tex",sep="")
	filenamecsv <- paste(tabDir,"table4.csv",sep="")
	latex(Table3,file=filename,rowname=NULL,caption=cap, where="!htbp")
  cat(paste("Saved table ",filename,"...\n",sep=""))
  cat(paste("Saved table ",filenamecsv,"...\n",sep=""))
}

table.5 <- function(){	#proportion at age data used in the assessment.
	aa <- round(A$Adata[,1:14]*100,2); colnames(aa)=paste(2:15)
	x <- data.frame(Year=A$ayr,age=aa) 
	cap <- "Age-composition (reported in percentages) of the combined U.S. and Can. 
		      commercial catch from 1977-2009. Age-15 represents a plus group."
	filename <- paste(tabDir,"Table4.tex",sep="")
	latex(x,file=filename,rowname=NULL,caption=cap,landscape=T, label="table.Adata", where="!htbp")
	
	aa <- round(A$page*100,2); colnames(aa)=paste(2:15)  
	x <- data.frame(Year=A$iyr, age=aa)  
	cap <- "Age-composition (percent) from acoustic surveys from 1977-2009.  
		      Note that age-15 represents a plus group. Proportions-at-age were
		      constructed by multiplying the conditional age-length data by the
		      length frequencies and collapsing over each size interval."      
	filename <- paste(tabDir,"Table5.tex",sep="")
	latex(x,file=filename,rowname=NULL,caption=cap,landscape=T,label="table.Page",  where="!htbp")
	print("Updated table 5")
}

table.cwa <- function(){	#sideways table for the catch-weight at age
	cap <- "Assumed mean weights-at-age in the commercial catch."
	xx <- round(A$cwa,3)
	colnames(xx) <- c("Year",paste("age",2:15))
	filename <- paste(tabDir,"CWA.tex",sep="")
	latex(xx,file=filename,rowname=NULL,caption=cap, size="footnotesize",label="TableCWA")
	print("Updated table cwa")
}

table.forecast <- function(){
  fn <- paste(admbDir,"tinss.for",sep="")
	dd <- read.table(file=fn, head=T)
	iyr <- unique(dd$Year)
	ctv <- unique(dd$CtStream)
	dtable <- NULL 
	quan <- c(0.05, .25, .5, 0.75, 0.95)
	
	for(j in ctv){
		d <- subset(dd, CtStream==j)
		for(i in iyr){
			dep <- quantile(subset(d, Year==i)$depletion, na.rm=T,quan)
			sbt <- quantile(subset(d, Year==i)$Sbt, na.rm=T, quan)/2
			fspr <- quantile(subset(d, Year==i)$fspr, na.rm=T, quan)
			f40spr <- quantile(subset(d, Year==i)$f40spr, na.rm=T, quan)
			ABC <- unique(subset(d, Year==i)$ABC)
			iABC <- quantile(subset(d, Year==i)$ABC, na.rm=T, 0.5)
   			OY <- quantile(subset(d, Year==i)$OY, na.rm=T, 0.5) 
 		  #print(round(c(i, fspr), 2))
			dtable <- rbind(dtable, c(i,iABC*1e6,OY*1e6,round(sbt, 2),round(dep, 4),round(fspr, 2), round(f40spr, 2)))
			#dtable=rbind(dtable, c(i,iABC*1e6,OY*1e6,round(sbt, 2),round(dep, 4),round(fspr, 2), round(f40spr, 2)))
		}
	}
	#US decision table
	ic <- c(1:13, 19:23)	
	filename <- paste(tabDir,"USDecisionTable.tex",sep="")
	cap <- "Decision table with three year projections"
	colnames(dtable) <- c("Year","ABC","OY",rep(paste(c(5, 25, 50, 75, 95), "th", sep=""), 4)) 
	dum <- latex(dtable[, ic],file=filename, rowname=NULL, caption=cap, label="tableUSdecision")
	write.table(dtable[, ic], file="USDecisionTable.csv", sep=",", row.names=F)
		
	#CAN decision table
	ic <- c(1:3, 14:18, 19:23)	
	filename <- paste(tabDir,"FDecisionTable.tex",sep="")
	cap <- "Decision table with three year projections"
	colnames(dtable) <- c("Year","ABC","OY",rep(paste(c(5, 25, 50, 75, 95), "th", sep=""), 4))
	write.table(dtable[, ic], file="CanDecisionTable.csv", sep=",", row.names=F)
		
	dum <- latex(dtable[, ic],file=filename, rowname=NULL, caption=cap, label="tableFdecision")
	 
	print(dtable)
	return(dtable)
}

table.mle.fore <- function(){
  fn<-paste(admbDir,"tinss.forecast",sep="")
	dd <- read.table(fn,header=T)
	dd[, 2:3] <- dd[, 2:3]*1e6
	dd[, 4:7] <- round(dd[, 4:7], 2)
	filename <- paste(tabDir,"ForecastTable.tex",sep="")
	cap <- "Three year projections of maximum likelihood-based Pacific Hake ABC,  OY, 
	        female spawning biomass,  spawning biomass  depletion level,  and relative SPR values
	        based on the 40:10 harvest control rule with F40 (top three rows) and fmsy 
	        (bottom three rows) overfishing targets."
	print(dd)
	print("Updated table mle.fore")
}

mb.table.all.param.est <- function(roundDec=2,formatOut="%1.2f"){
  mcmc <- A$mc
  mcmc <- mcmc[(Burn+1):nrow(mcmc),]

  numParams <- 2
  d <- as.data.frame(matrix(nrow=numParams,ncol=2))
  colnames(d) <- c("Parameter","MCMC median")
  row <- 1

  d[row,1] <- "Log(R0)"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$log.ro),roundDec))
  row <- row + 1
   
  d[row,1] <- "h"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$h),roundDec))
  row <- row + 1
 
  d[row,1] <- "Log(M)"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$m),roundDec))
  row <- row + 1

  d[row,1] <- "Log of average recruitment"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$log.rinit),roundDec))
  row <- row + 1

  d[row,1] <- "Variance ratio (rho)"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$rho),roundDec))
  row <- row + 1

  d[row,1] <- "Inverse total stdev (varphi)"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$varphi),roundDec))
  row <- row + 1

  d[row,1] <- "Survey Age at 50% 1st harvest"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$ahat.surv),roundDec))
  row <- row + 1

  d[row,1] <- "Survey Age at 50% 1st harvest stdev"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$ghat.surv),roundDec))
  row <- row + 1

  d[row,1] <- "Commercial Age at 50% 1st harvest"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$ahat.comm),roundDec))
  row <- row + 1

  d[row,1] <- "Commercial Age at 50% 1st harvest stdev"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$ghat.comm),roundDec))
  row <- row + 1

  d[row,1] <- "Survey catchability (q)"
  d[row,2] <- sprintf(formatOut,round(median(mcmc$q),roundDec))
  row <- row + 1

  years <- ryr

  # add F parameters
  burnedF <- log(A$mc.ft[(Burn+1):nrow(A$mc.ft),])
  tmpMeds <- apply(burnedF,2,median)
  currYear <- yrs[1]  # starts with 1966 for Hake
  for(year in 1:ncol(A$mc.ft)){
    d[row,1] <- paste("Log fishing mortality #",year," - ",currYear,sep="")
    d[row,2] <- sprintf(formatOut,round(tmpMeds[year],roundDec))
    currYear <- currYear + 1
    row <- row + 1
  }   

  
  # A$mcRecDevs should be 60 columns long for 2012, 1952-2012.
  burnedRec <- A$mcRecDevs[(Burn+1):nrow(A$mcRecDevs),]
  tmpMeds <- apply(burnedRec,2,median)
  currYear <- years[1]-nage  # starts with 1952 for Hake
  for(year in 1:ncol(A$mcRecDevs)){
    d[row,1] <- paste("Mean recruitment #",year," - ",currYear,sep="")
    d[row,2] <- sprintf(formatOut,round(tmpMeds[year],roundDec))
    currYear <- currYear + 1
    row <- row + 1
  }   

  fn <- paste(tabDir,"mb.table.est.parameters.csv",sep="")
  write.csv(d,file=fn,row.names=F)    
}
